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ABSTRACT 

In many strategic shallow water areas the geoacoustic properties of the 
sub-bottom are largely unknown. In this thesis it is demonstrated that 
inverse theory and measured data from a single hydrophone can be used to 
accurately deduce the geoacoustic properties of the sub-bottom, even when 
the initial background geoacoustic model is a highly inaccurate “guess”. 
Since propagation in shallow water is very sensitive to the geoacoustic 
properties of the sub-bottom, the inverse technique developed in this thesis 
presents the Navy with a vitally important, practical, and inexpensive means 
to improve sonar performance prediction in a potentially hostile environment. 

To provide ground truth for the inverse technique, measured data 
collected during Project GEMINI were compared to the inverse solutions. 
Detailed, site-specific geoacoustic models were developed for two array 
locations and the Finite Element Parabolic Equation (FEPE) model was used 
to estimate transmission loss (TL). The model estimates from FEPE 
compared well with the measured data and the detailed geoacoustic models 
were considered as “ground truth”. To test the efficacy of the technique, 
initial background geoacoustic models were constructed assuming no a priori 
information of the bottom. The resultant inverse solution was used to predict 
the geoacoustic properties at each of the sites. The final results were in 
excellent agreement with the measured data and the resulting inverse 
technique TL estimates were as good or better than the TL estimates 
obtained from the detailed, site-specific geoacoustic models. 
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[. INTRODUCTION 


A. NEW PRIORITIES 

The United States Navy (USN) strategic vision articulated in "...From 
the Sea,’’(O'Keefe et al., 1992) asserts U.S. naval forces as full participants in 
the new national strategy for the post-Cold War Era. No longer focused 
primarily on open ocean global conflict, the naval service’s new emphasis is 
clearly placed on regional contingencies in the littoral regions of the world. 
Accordingly, the Navy has downgraded its previous anti-submarine warfare 
(ASW) requirements that stressed blue-water capabilities against Soviet 
attack and ballistic missile submarines to a more localized focus towards 
conventionally powered submarines operated by Third World countries. 

Apart from the United States and the Commonwealth of Independent 
States (CIS), over 40 other countries now operate more than 400 submarines 
worldwide (Withers, 1992), half of which are operated by Third World 
countries. Some of the more modern submarine types held in Third World 
inventories include: the Ex-Soviet KILO-Class, the German Type 209, the 
British OBERON-Class and the French AGOSTA-Class. Historically, diesel 
electric submarine (SSK) operations have been limited in both endurance and 
range. However, recent developments in propulsion technology, including air 
independent propulsion (AIP), have substantially reduced snorkeling 
requirements and have also produced a quieter, more formidable SSK. Due 
to their smaller size, SSKs offer lower target strengths than their nuclear 
counterparts; they offer less of a target to ping on and thus produce less 
return. When operating slowly, they are non-cavitating and provide little 
Doppler return. The net effect of these improvements is to reduce the 
detection opportunities for ASW forces. In order to gain control of the seas, 
whether denying it to an adversary or protecting it for self use, the U. S. 
Navy must develop the ability to conduct successful ASW operations within 
the littoral environments. 



B. BACKGROUND 

Shallow water acoustics has long been a topic of ASW research and 
has been thoroughly investigated both theoretically and experimentally (e.g., 
Akal, 1980; Frisk et al., 1984; Rajan et al., 1987). However, the accumulation 
of theory and direct measurements has failed to provide consistent, accurate, 
and qualitative prediction of acoustic propagation in shallow water 
environments. The reason for this is due almost entirely to the complexity of 
the problem. In shallow water the acoustic medium consists of properties 
that vary both spatially and temporally on rather short scales. Since the 
spatial variation of geoacoustic properties of the bottom/sub-bottom is large, 
it is unlikely that a complete and accurate data base of bottom/sub-bottom 
geoacoustic properties will ever be collected for all shallow water areas of 
interest. Therefore, the U. S. Navy critically needs and accurate, cost 
effective, and practical method to measure or deduce the geoacoustic 
properties of the sub-bottom. The aim of this thesis is to solve this 
challenging problem vital to the U. S. Navy’s ability to accurately model 
acoustic propagation in shallow water. 

Unlike the deep-water environment, acoustic propagation in shallow 
water is dominated by repeated interactions with the seafloorand sub-bottom. 
A distinguishing characteristic of shallow water propagation is that the 
sound-speed profile (SSP) is generally nearly constant with depth or is 
downward refracting, meaning that long-range propagation takes place 
exclusively via bottom-interacting paths. The deep-water concept of direct 
path (DP) and convergence zone (CZ) acoustic propagation does not apply in 
the shallow water environment. Instead, the more important propagation 
paths are of the form of refracted bottom-reflected or surface-reflected 
bottom-reflected propagation. Generally, the high attenuation resulting from 
repeated acoustic interaction with the bottom and limited water depths 
severely impacts propagation, thereby reducing detection ranges in shallow 
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water. Figure 1 gives an example of the degree of transmission loss (TL) 
variability as a function of range in shallow water. The plot illustrates 
average TL versus range for shallow regions (100-200 m) around the world 
(Urick, 1979). Two features are of immediate interest in this figure. One is 
the spread of the TL data in excess of 50 dB at 100 km due primarily to the 
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Figure 1. Transmission loss variability in shallow water (from Urick, 1979). 


varying geoacoustic makeup of the seafloor. The second feature is that 
propagation is generally better than that experienced with spherical 
spreading (20 log r) at short and intermediate ranges but is larger for 
extended ranges. This feature is attributed to reduced spreading loss due to 
the trapping of acoustic energy in the shallow water wave guide which tends 
to decrease TL at shorter ranges. However, due to increased interaction with 
the boundaries, TL increases at longer ranges (Jensen et al.. 1994). The 
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cross-over range depends on water depth, acoustic frequency, and the nature 
of the sea bottom. 

At low frequencies acoustic propagation in shallow water is highly 
dependent on the geoacoustic properties of the sea bed. As sound propagates 
through the sea bed, the acoustic pressure field is influenced by the unique 
physical properties of the sedimentary sequence through which it travels and 
also by the nature of the lithologic boundaries. Thus, an area of considerable 
interest in underwater acoustics is the determination and qualitative 
description of these geoacoustic parameters that effect acoustic propagation. 
The traditional approach for describing the acoustic interaction with the 
bottom has been to develop a geoacoustic model (direct method) by methods 
described in Hamilton (1980). Provided adequate geologic and geophysical 
data exist in the study area, this method has proved invaluable to acoustic 
propagation modeling efforts. 

A second approach involves the inversion of these properties from a 
measured sound pressure field. Numerous techniques have been cited in the 
literature for obtaining inverse geoacoustic solutions. Frisk et al. (1986) used 
amplitude versus range data obtained using a deep-towed pulsed continuous 
wave (CW) source and two receivers anchored near the bottom to infer 
geoacoustic properties of the upper sedimentary layers in deep water. The 
inferred geoacoustic model was derived from the measured data using an 
iterative forward modeling technique. The method relies on time-gating to 
separate the bottom reflected arrivals from the surface reflected signal. Zhou 
and Zhang (1987) deduced geoacoustic properties through normal mode 
filtering techniques. Their technique utilizes dispersion analysis and normal 
mode measurements to determine acoustic attenuation and propagation 
speeds for a horizontally stratified bottom as a function of frequency. This 
type of experiment requires the use of a vertical array of hydrophones to 
successfully separate the modal arrivals. 



In this thesis we use a perturbative inversion technique (Lynch et al., 
1987; Rajan, 1992) to obtain the compressional sound speed as a function of 
depth for the upper sedimentary layers using as input data the eigenvalues 
of the propagating sound field. The method consists of measuring the 
magnitude and phase of the pressure field as a function of range generated 
by a CW point source and numerically Hankel transforming these data to 
obtain the depth-dependent Green’s function versus horizontal wavenumber. 
The resulting Green's function contains information about the discrete and 
continuous modal spectra of the waveguide; specifically the prominent peak 
in the Green's function occurs at wave numbers corresponding to the 
eigenvalues for the trapped modes of the waveguide. Although the 
sedimentary layers have enough rigidity to support shear wave propagation. 
Fryer (1978) indicates that the effect of compressional to shear wave 
conversion at the layer interfaces is minimal for frequencies above 20 Hz. 
Thus in this work, the sediments are treated as a fluid extension of the water 
column and the shear wave effects are neglected. Furthermore, the water 
column is characterized as isospeed bounded at the surface and bottom by a 
horizontally stratified media. 

C. PROJECT GEMINI 

A series of narrow-band bottom interaction experiments, collectively 
called Project GEMINI, were performed in the Gulf of Mexico near Corpus 
Christi, Texas in September 1985. Project GEMINI was designed to be a 
“baseline" experiment to evaluate low-frequency acoustic interaction with 
the seabed (Lynch et al., 1991). The intent of the experiment was to obtain 
high-resolution acoustic measurements at a relatively benign (i.e., 
geologically simple, range-independent) acoustic environment which could 
later be used to test or benchmark state-of-the-art acoustic propagation 
models. In order to properly address the performance and robustness of 
propagation models in shallow water, the models must first be evaluated in 
5 



regions where the physical parameters (e.g., sound speed, geologic structure, 
etc.) affecting acoustic propagation are well understood. A propagation model 
that does not perform well in a geologically simple area is less likely to 
perform well in more geologically complex regions. The GEMINI region is 
characterized by a smooth, gently sloping seafloor overlying a sequence of 
horizontally-layered sediments. Physical characteristics of the GEMINI 
seafloor such as sediment thickness, speed, and density are well-known and 
documented (Matthews et al., 1985). Thus, the experiment provided the 
opportunity to qualitatively determine how well one could predict measured 
acoustic data using propagation models with high quality geoacoustic 
information as input to the models. 

Project Gemini consisted of five synthetic aperture array experiments 
occupying three different array locations along the continental shelf of 
Texas. The experiments consisted of towing a narrow-band CW source at a 
fixed depth with output tonals of 50 Hz and 140 Hz, away from a pair of 
moored receivers (Lynch et al., 1991). The track lines in each of the 
experiments extended from zero range at the array locations to a maximum 
distance of about 5000 m. The lines were oriented such that they transversed 
upslope from the array positions nearly perpendicular to the isobath lines. 
The two moored receivers were placed midway in the water column and 
approximately 1.5 m above the bottom, respectively. A complete listing of 
array locations, source and receiver combinations, and dates of the 
experiments are presented in Table 1. 

The primary array site selected for Project GEMINI, here referred to 
as the Rubano Site, was previously studied by the Applied Research 
Laboratory, Pennsylvania State University (ARL/PSU) in 1975 (Rubano, 
1980). Two alternate sites, a shallow-site to the north and a deep-site to the 
southeast were also incorporated into the exercise (Figure 2). The water 
depth for the shallow site is approximately one-half that of the Rubano site 

6 



whereas the deep-site is about twice the depth. Prior to the exercise a 
detailed environmental assessment of the GEMINI region was carried out by 
researchers at the Naval Research Laboratory-Stennis Space Center 
Mississippi (NRL-SSC) (Matthews et al., 1985). This environmental 
assessment produced an acoustical evaluation of the test site which used as 
input a high-resolution geoacoustic model. The geoacoustic model was 
constructed from the most current multi-channel seismic, bathymetric, and 
drill-hole information available at the time. The geoacoustic model was 
specifically tailored for the Rubano site although it was hoped, given the 
proximity of the array locations and the presumed benign geologic properties 
of the region, that the model could be extrapolated to the two other sites. 


Experiment 

Number 

Date 

Water 

Depth <m) 

Source Depth 

(ml 

Receiver 

Depths (m) 

Description 
































Table 1. Description of Project GEMINI. 










Figure 2. Geographic locations of the Project GEMINI experiments off the Texas 
Coast. Isobaths in meters. 

D. OBJECTIVE 

In shallow water the geologic and geoacoustic properties of the bottom 
can have significant variations over rather small ranges. Current 
geoacoustic modeling methods (point models) only address the physical 
properties of the sediment and rock as a function of depth below the sea floor 
at a single point and not as a function of horizontal distance. In order to 
adequately account for the spatial variations in geologic characteristics a 





series of geoacoustic models must be developed to account for this variability. 
However, as stated previously, such spatially varying data are not generally 
available. An alternative approach is to infer these geoacoustic parameters 
through geoacoustic inversion techniques. 

The objectives of this study are three-fold. The first objective is to 
compare several propagation model estimates of TL with measured low- 
frequency acoustic propagation data collected during Project GEMINI. Such 
a comparison will determine which of several TL models most closely 
approximates the observed TL pattern, i.e., which model solves the forward 
problem best. The acoustic models will have as input location-specific, high- 
resolution geoacoustic characterizations of the seafloor. The acoustic models 
that will be used in this study are the SACLANTCEN normal mode model 
(SNAP), the NRL normal mode model (KRAKEN), and the NRL finite 
element parabolic equation model (FEPE). 

The second objective is to infer the geoacoustic properties of the 
seafloor in the GEMINI region using the perturbative inversion method. The 
underlying assumption of perturbative inversion theory is that a reasonable 
"first guess'' of the unknown geoacoustic properties must be incorporated as a 
starting point in the computation process (Frisk et al., 1989). L’sing a 
detailed, location-specific geoacoustic model as a starting point in the 
inversion process, it will be shown that a geoacoustic model, derived through 
inversion, results in improved TL model estimates over those presented 
earlier. 

Lastly, an objective vitally important to the Navy’s operational 
performance in shallow water will be pursued. In many shallow water 
regions little geologic or geophysical information is readily available to 
formulate useful geoacoustic characterizations of the bottom. If. however an 
“educated first guess" for a geoacoustic model is made in a shallow water 
area of interest, it is important to determine how accurate the perturbative 
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inversion technique will be in estimating the actual geoacoustic properties. 
Moreover, the TL model results using this process must be assessed to see 
how comparable they are to those formulated using a high quality 
geoacoustic model as an input. 
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II. NORMAL MODE SOLUTIONS IN A WAVEGUIDE 


A. MODAL PROPAGATION 


The inverse method considered in this thesis requires as input data 
the measured horizontal wave numbers or eigenvalues which characterize 
the modal propagation of sound through the ocean and seafloor. Since the 
ideas presented in normal mode theory are essential to the inversion process, 
a brief discussion of normal mode propagation is presented. Normal mode 
theory provides a ‘full wave’ solution of the wave equation in the ocean 
waveguide which takes into account the physical characteristics of the 
medium, frequency dispersion effects, and interaction of the sound field at 
the boundaries (Casey, 1988). 

The propagation of monochromatic sound in the ocean is described 
mathematically by solutions to the reduced wave or Helmholtz equation. 
Using the assumption of cylindrical symmetry, the spatial part of the 
acoustic pressure p(r,z,Zo), due to a point source at range r = 0 and depth z = 
Zo with a sinusoidal or harmonic time dependence exp(-icot), satisfies the 
Helmholtz equation 


|_r drVdrJ dz 2 


, , 1 -28(z-zJ 

k 2 (z) p(r,z,z 0 ) =-8(r) 


( 1 ) 


where p(r,z,z<,) represents the acoustic pressure at range r and depth z, co is 
the radian frequency, and k(z) = co/c(z) is the total wave number. 

Appropriate boundary conditions must be specified in order to 
determine a unique solution to Equation (1). At the interfaces where there 
are significant variations in properties, e.g., at the water sediment interface, 
continuity of pressure and normal particle velocity must be met. The 
symmetry and range independence of Equation (1) leads to the method of 
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separation of variables in solving the equation. By expressing the pressure 
field p(r;z,Zo) by 


p(r,z,Zo) =Z(z)R(r) (2) 

and substituting this into Equation (1) leads to the following expression 

ifi d f dnY! if d 2 z] 
rLt drv dr Jj + z[_*^ Z dz 2 l H 


+ k 2 = 0 


(3) 


The first term in Equation (3) is a function of r only whereas the 
second term is a function of z only. This equality is satisfied if and only if 
both expressions are equal to the same constant of proportionality which 
leads to the two separated equations: 


d 2 R(r) 1 dR( r) 
dr 2 r dr 


k 2 R(r) = 0 


(4) 


and —k 2 Z(z) = 0 (5) 

dz 

where (k r ) 2 and (k,) 2 are the separation constants defined by 

k 2 = k 2 (z)-k 2 (6) 


Equation (4) is a function of range only, and thus from Equation (6) k r , the 


horizontal component of the wave number, must be constant throughout the 
waveguide. This can also be expressed in terms of Snell’s law since 


k = ksin0 = — sin0 (7) 

c 

where 0 is the vertical angle of incidence and co is constant, and (sin9)/c is 
constant (Clay and Medwin, 1977). Accordingly, at a zero incidence angle 
(0 = 0°), there is no horizontal contribution to the total wavenumber. The 
solution to Equation (4) is a cylindrical Bessel or Hankel function of the 
zeroth order. 

R(r) = H<"(k r r) (8) 

For distances of more than a few wavelengths from the source, the 
asymptotic expansion of H^’fkjr) (Haberman, 1987) can be expressed as 
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The exponential expansion of the cosine term in Equation (9) yields the 
following characteristic equation 

> J^'4 ik " r -f} +esF (r ,k - r+ f]j <10) 

When associated with a harmonic time dependent term, Equation (10) 
represents the total wave motion consisting of an inward propagating wave 
represented by the first term on the right side of Equation (10) and an 
outward propagating wave represented by the second term on the right of 
Equation (10). Since we are only interested in an outward propagating, 
exponentially decaying waveform, only the second term is used in the 
formulation (Clay and Medwiri, 1977). 

lu> 

1. Rigid Bottom Solution 

Assuming that the ocean surface represents a pressure release 
surface and that the ocean bottom at z = h is perfectly rigid, allowing no 
acoustic penetration, the “normal mode” differential equation (Equation (5)) 
and the accompanying boundary conditions are: 


^ + k:z,„=o 

dz 

(12) 

dZ(h) 

Z(0)=0, = 0 

dz 

(13) 

By solving this equation as a Sturm-Liouville eigenvalue 

general solution to Equation (12) is of the form 

problem, the 

Z(z)=Cisin(k*z)+C2COs(k z z) 

(14) 

where Ci and Cs are arbitrary constants and k; is the vertical wavenumber. 

The boundary condition at the free surface (z=0) requires that the acoustic 
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pressure to be zero for all range, r, and time, t, thus requiring C2=0. 
Furthermore, an incident plane wave on the sea surface is assumed to be 
completely reflected, undergoing a 180° phase shift upon interaction with the 
free surface. The rigid bottom boundary condition is characterized by the 
constraint that the acoustic pressure is a maxima at the seafloor at z=h and 
that the normal particle velocity vanishes on the surface resulting in 
complete reflection (Casey, 1988). This means that k z z = n/2, 3 tt/ 2, 5it/2... or 
more generally 

(kz)m = (7i/h)(2m-l), m = 1, 2, 3 ... (15) 

(k 2 ), n is termed an eigenvalue and the solution is an eigenfunction (normal 
mode) of the problem. There exists an infinite number of eigenvalues and 
thereby an infinite number of eigenfunctions Z m (z) as solutions to this 
boundary value problem. Each solution is characterized by a mode shape 
Zm(z) and a vertical propagation constant (kz) m . The subscript m designates 
the modal number for a particular eigenvalue or eigenfunction. 

The pressure field in Equation (2) can now be written as 

p(r,2)=£E.(r)Z„(z) (16) 

Substituting this expression for pressure into Equation (1) yields: 


The eigenfunctions of the Sturm-Liouville problem form an orthonormal set 
of square integrable functions on 0<h<z such that 


r Z.(z)Z,(z) fO m * n 

j p(z) jl m = n 


(18) 


Multiplying both sides of Equation (17) by 



(19) 


and applying the orthogonality relationship of (18) yields the following 
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(20) 


The solution to Equation (20) is a Hankel function of the first kind such that 


R «W = T7-T z ra ( z -) H i"( k r *> (2D 

■lp(z„) 


Substituting the expression for H^’(k r r)(k n r) found in Equation (11), the 
solution to Equation (1) can be expressed as a sum of normal modes 


exp(ik r r) 



( 22 ) 


K 


Equation (22) describes the pressure field as a linear superposition of 
traveling cylindrical waves propagating outward from a source. Each of 
these modes is a standing wave in the vertical direction with a unique depth 
distribution described by the eigenfunction Z m (z) and the vertical modal 
wavenumber (k:)m The vertical wavenumber must ensure that a pressure 
node exists at the surface and, in the case of a perfectly rigid bottom, a 
pressure anti-node at the seafloor. The eigenfunctions Z m (z) essentially 
determine the depth dependence of the acoustic pressure field. For a given 
water depth and source frequency, only a finite number of modes may exist. 
For large mode numbers m. the horizontal wavenumber becomes imaginary 
and the solution defined by Equation (22) rapidly decays with increased 
horizontal ranges. 

2. Penetrable Bottom Solution 

If the ocean bottom is not a perfect reflector, sound energy will 
penetrate into the sea bed and solutions to the Helmholtz equation must be 
found both in the water column and the ocean bottom. The solution to 
Equation (1) for the penetrable bottom case can be obtained by utilizing the 
Green's function technique. The conditions imposed upon this case are the 
same as that of the rigid bottom with the exception that the velocity 
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potential is discontinuous at the seafloor. A brief derivation of the Green’s 
function is now given. 

Assuming a constant density water column of thickness (h), sound 
speed c(z), bounded at the surface and bottom by a horizontally stratified 
medium (Figure 3), the spatial part of the sound pressure field (p) due to a 
harmonic time dependent point source located at r = 0 and z = Zo satisfies 
Equation (1). 

The solution to Equation (1) can now be written as the zero-order 
Hankel transform of the depth-dependent Green’s function g(k t ,z). This is 
achieved by taking the zero-order Hankel transform of both sides of Equation 
(1) where the Hankel transform pair is defined by the following: 


p(r,z,) = Jg(k,,z)J 0 (k [ r)k,dk, 

(23) 

g(k r ,z) = J p(r)J # (k,r)rdr 

(24) 


The Hankel transform variable, k r , transforms the solution from a 
depth and range space p(z,r) to a solution as a function of depth and 
horizontal wave number space (z, k r ) (Evans, 1975). 
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Figure 3. Horizontally stratified ocean model. 




After applying Equation (24) to Equation (1) and re-ordering terms, 


Equation (1) is transformed into the depth separated wave equation 

|"^,+k 2 (z)-k r ] gk . z z ) = -25(z- z 0 ) (25) 

dr r 0 


Here g(k r ,z) is the depth-dependent Green’s function, J u is the zero-order 
Bessel function, r is horizontal range, and k, is the horizontal wave number. 
The general solution to Equation (25) can be written as the sum of the 
particular solution (non-homogeneous part) and the characteristic solution 
(homogenous part). The characteristic solution of Equation (25) is of the form 
g(k r ) = <|>s.<t>B (26) 

where <|>k and dn are linearly independent solutions of Equation (25) chosen 


such that te satisfies the surface boundary condition, and <|>b satisfies the 
bottom boundary condition. By applying the method of variation of 
parameters, the Green's function can be defined by: 


g(kr.z) 

g(k r ,z) 


-2iji s (k r ,z)d n (k r ,z 1 ,) 

W(zJ 

- -2^(k,,z „ )de(k,,z) 

W(z 0 ) 


0 <z< z 0 
z t Sz5h 


(27) 

(28) 


such that W(Zo) is the Wronskian given by 

W(Zo) = <1 »b(Zo) l)l's(Zo)- d' b(Zo) <t>s(Zo) (29) 

where the prime indicates differentiation with respect to z. <J>s and ()>b also 
must satisfy the following impedance relations. 



=n and are input quantities that contain boundary conditions necessary for 
the calculation of the Green’s function (Frisk et al., 1980). Both quantities 
are dependent on the horizontal wavenumber, frequency, and the acoustic 
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properties of the bottom and the surface. 4 b depends on the water density, 
p«. and the sound speed c(z), at the water-bottom interface, while 4s depends 
on pxv and the sound speed c(z) at the water-surface interface. Equation (30) 
therefore can be expressed in terms of plane wave reflection coefficients, Rq 
and Rs. for the half-space corresponding to sound speed cb and cs bounding 
the bottom and the surface, respectively. 

<j>B(z) = Ci[exp(-ik z z)+ RB(kr) exp(ik z z)] (31) 

<J>s(z) = C2[exp(-ik z z)+ R,(k r ) exp(ik z z)] (32) 


where k z 2 +k r 2 = k 2 =(to/c) 2 is the total wavenumber and Ci, C2 being arbitrary 
constants. Equation (31) can also be re-written as 



i(l + R B ) 
k z (l - R b ) 


at z - 0 


_<t> i _ -i l + R s exp(-2ik z h) 

<t>' % k 2 |_l-R s exp(-2ik z h) ’ 


(33) 


(34) 


Substituting these expressions into Equations (27) and (28) and 
further assuming that the ocean surface is a pressure release boundary (i.e., 
Rs = -1), the following expression for the depth dependent Green’s function 


can be found: 
g(kr,z) = 

where z+ - z+Zo, 


+ R B (k [ )e 2 ‘ k ' 


-ik z Jl + R B (k [ )e 2,k 


In order to evaluate the integral expression for p(r,z) the Bessel 
function, J a (k f r) in Equation (23), is expressed in terms of Hankel functions 



J„(k,r) = ^-[Hf, ll (k r r) + H!, : ’(k r r)] 


(36) 


The following relation 

[HI,’’ e.\p(-ijtk r r) = -H' 1, (k I r)l (37) 

allows the range of integration of the transform integral to be evaluated over 
the entire real k r axis from -oo to co (Casey. 1988). The solution of Equation 
(23) is then found through contour integration methods (Evans, 1975) and is 
expressed as the sum of discrete and continuous portions of the normal mode 
spectrum given by 

p(r,z).m^Z,(z«.(z.)Hi"<k,_r)+IW (38) 

where the eigenvalues k ra and eigenfunctions Z m satisfy 

[£ r + k , (z)-k;_}z,„(z)=0 (39) 

The discrete sum in Equation (38) corresponds to modes perfectly 
trapped in the waveguide and is of the same form as that of the rigid bottom 
solution. Typically, this trapped mode sum dominates the long range 
behavior of the acoustic field whereas the continuum contribution I(r) is only 
significant at short ranges (Lynch et al., 1987). The continuum is rapidly 
attenuated with increasing range from the source and is typically neglected 
for ranges over a few wavelengths from the source. The eigenvalues for the 
penetrable bottom case have the same interpretation as those in the rigid 
bottom solution. However, they will take on different values since sound is 
propagating in both the water column and the sub-bottom and the 
characteristic equation will be different. 

The modal representation is closely correlated with the analytic 
properties of the Green’s function. A knowledge of the behavior of g(k r ) (for 
example, the number of resonances and their spatial location in regard to 
horizontal wavenumber) provides information about the nature of the modal 
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propagation in the waveguide. The Green’s function is characterized by a 
finite valued continuum with a discrete set of resonances occurring when the 
denominator of Equation (34) is set equal to zero. These poles of g(k r ,z) 
correspond to the eigenvalues of the trapped and virtual modes excited in the 
waveguide. The position and magnitudes are highly influenced by the 
inherent acoustic properties of the bottom (Frisk et al., 1986) and the 
geometry of the waveguide. As a result, the basic principle underlying the 
inverse method presented in this paper is first to obtain an estimate of g(k r ) 
from measurements of p(r) by numerically performing the Hankel transform 
in Equation (23). The modal features of g(k r ) are then used to determine a 
geoacoustic model. A specific perturbation technique relating the 
eigenvalues of the trapped modes to the geoacoustic model are given in the 
next section. 

B. PERTURBATIVE INVERSION TECHNIQUE 

A detailed discussion of the perturbative inversion theory that is the 
basis for the analysis in this thesis is presented in Rajan et al. (1987). The 
objective of the perturbative inversion technique is to estimate the variation 
in geoacoustic parameters (i.e., sound speed, density, attenuation, etc.) from 
the differences between the measured eigenvalues and the eigenvalues 
computed for a background geoacoustic model (Frisk et al., 1986). In the 
perturbative technique, an initial background geoacoustic model is assumed 
and the differences of the geoacoustic parameters from the geoacoustic 
parameters in the initial background model are calculated. For the purpose 
of this thesis, we are concerned with determining only the compressional 
sound speed profile in the bottom, c(z), for the deep and shallow GEMINI 
sites. Perturbative techniques to determine other geoacoustic parameters 
(i.e., density, attenuation, etc.) are not addressed. However, these 
parameters can be modified to be geologically consistent with the final 
estimate of c(z). As a first step, an initial background geoacoustic model for 
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c(z) is assumed, and a major objective of this thesis is to evaluate the 
perturbative technique when the initial estimate of c(z) must be based on 
little to no geoacoustic information about the area. Perturbation theory then 
provides the means for determining the difference Ac(z) between the 
background model and the true earth model. 

From first-order perturbation theory, an integral equation relating the 
perturbation of eigenvalues is given by: 



(40) 


where Akm represents the difference between the experimentally measured 
modal eigenvalues and the eigenvalues computed from the background 
geoacoustic model, km and Z m (z) are the with eigenvalues and normalized 
eigenfunctions for some assumed background geoacoustic model, p(z) is the 
density profile for the background geoacoustic model, k(z) is the acoustic 
wavenumber defined by co/c(z) (where co is the angular frequency and c(z) is 
sound speed profile of the geoacoustic mode;), and Ac(z) is the difference in 
sound speed between the background and true models. All terms in 
Equation (40) are known except for Ac(z) which is the solution being sought. 

Equation (40) is a Fredholm integral equation of the first kind for Ac(z) 
and generally does not produce a unique solution unless constraints are 
made upon its solution (Lynch et al., 1987). Reordering the known quantities 
in Equation (40) such that 



(41) 


results in 



(42) 
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K(m,z) is defined as the kernel of the integral equation and consists of the 
known physical parameters for a discrete mode and depth step. Using 
methods proposed by Backus and Gilbert (1970), a solution for Ac(z) is given 
by: 

Ac(z) = XXK(m,z)Ak m A"'(m,n) (43) 

where A(m,n) is the spectral decomposition of the kernel K(m,z) and 

A(m,n) = }K(m,z)K T (n,z)dz (44) 

The solution given by Equation (43) is called the Backus-Gilbert 
solution for perfect data. This solution represents the most “Dirac delta like" 
resolution kernel and was used in this analysis because smoothing the 
solution was not desirable until after c(z) was estimated. This represents a 
solution without constraints and is considered a fundamental property of the 
perturbative technique not addressed originally in Rajan et al., (1987). If 
constraints are used in the solution for c(z), the most natural (i.e., 
unconstrained) solution characteristics will not be realized. The “price to 
pay” for the natural or “unconstrained” approach is instability in the solution 
when the ratio of the highest to lowest eigenvalue is large. In that instance, 
it is well known that the pseudo-inverse to the matrix in Equation (44) is ill 
conditioned, but a stabilization technique based on the Frobenious norm was 
used. Also, a simple “box car” average of c(z) in depth was performed to 
dampen fluctuations in c(z) with depth. 

As a first step to the inversion process, measured complex pressure 
field versus range data was transformed from the range domain to the 
horizontal wave number domain (depth-dependent Green’s function) via the 
Hankel transform pair. As discussed earlier, the normal modes of the 
waveguide show up as strong resonance peaks. The measured modal 
eigenvalues are experimentally measured by simply picking the 
peaklocations in \g(k r ) 1 (Lynch et al., 1991). Figure 4 is a typical Green’s 
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function versus horizontal wavenumber plot for a particular environment. In 
this instance, each of the peaks (modes 1-4) has a characteristic wavenumber 
or eigenvalue identified in our formulation as kmtm 


Green's Function versus Horizontal Wavenumber 
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Figure 4. Plot of the depth-dependent Green’s function | g(kr) | versus horizontal 
wavenumber k at 140 Hz. The normal modes of the particular waveguide are 
identified by poles of g(kr). 


The NRL normal mode model, KRAKEN, was used to calculate the 
modal eigenvalues and eigenfunctions, km and Z m (z), respectively, using the 
assumed background geoacoustic model. In this section we assume a 
background geoacoustic model with a linear sound speed gradient and with 
constant density and attenuation with depth to illustrate the basic concepts 
of the inversion technique. In order to demonstrate the effectiveness of the 
perturbative inversion technique, the initial sound speed at the water-bottom 
interface and the sound speed gradient were also varied. Once km and Z m (z) 
are calculated, Ak m is found by the relation 

Ak m =k rae »,-k m (m=l,2,3...) (45) 

These values were entered into Equation (41) and the result, Ac(z ), is 
calculated by Equation (43) and represents the perturbation of the sound 
speed profile about that proposed in the initial background model. The sound 
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speed profile is then adjusted by subtracting Ac(z) from the initial value of 
c(z). The perturbative process continues until Ac(z) is minimized. 


1. Synthetic Test Cases 

In this section we will apply the perturbative inversion technique to 
a set of synthetic data for which the correct solution is known. The synthetic 
earth model along with the two initial background geoacoustic models are 
shown in Figure 5. The synthetic earth model consists of a 50 m thick 
isospeed water column (1500 m/s), overlying a 10 m thick isospeed sediment 
layer (1800 m/s), which comformably lies over a 60 m thick linear gradient 
(Vc =1.50 S' 1 ) sedimentary layer. The basement consists of an isospeed half¬ 
space with a compressional wave speed of 1900 m/s. The density is assumed 
to be constant p = 1.6 g/cm :i for all sediment layers and half-spaces. Although 
the water column sound speed profile for this example is assumed to be 
isospeed, it can consist of any structure. In general, the water column profile 
will be known from available XBT and CTD measurements and will thus 
drop out of the perturbation integrals, leaving only the sediment properties to 
be determined (Rajan et al., 1987). 
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The frequencies that will be considered in the synthetic examples are 
25, 50 and 150 Hz. Sub-bottom penetration by acoustic energy is highly 
dependent on frequency with lower frequencies being able to penetrate 
deeper into the seafloor than higher frequencies. As a result, lower 
frequencies are able to resolve the deeper sedimentary structure, whereas 
the higher frequencies resolve the finer near-surface structure. The synthetic 
Green’s function versus horizontal wavenumber plots for each of the 
frequencies are shown in Figures 6a-c These poles of g(kr,z), as discussed 
previously, correspond to the eigenvalues of the trapped and virtual modes 
that are excited for the particular waveguide at the given frequency. In 
general, higher frequencies excite more modes and thus have more trapped 
eigenvalues than do lower frequencies. The modal eigenvalues for the 
synthetic earth model (kmoa) are listed in Table 2. 

Two initial geoacoustic background models were assumed in this 
example: (1) an 1850 m/s isospeed model, and (2) a linear gradient model c(z) 
= 1750 + 2.0 s Kz). For the purposes of these test cases, density was taken to 
be a constant p = 1.6 g/cm 3 , and compressional attenuation was assumed 
negligible. The modeled eigenvalues (km) for the initial background models 
are shown in Tables 3 and 4. 
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Green’s Function versus Horizontal Wavenumber 
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Figure 6. Synthetic Green's function magnitude versus wavenumber for 25 Hz. 


Green's Function versus Horizontal Wavenumber 
(50 Hz) 


16000 
e 14000 

•f 12000 

| 3 10000 • 

^ | 8000 • 
a Jj 6000 
i — 4000 - 
° 2000 ■ 

0 J— 
0.17 


mode 2 





0.18 0.19 0.2 0.21 

Horizontal Wavenumber (1/m) 


Figure 7. Synthetic Green’s function magnitude versus wavenumber for 50 Hz. 
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Green's Function versus Horizontal Wavenumber 
(150 Hz) 
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's function magnitude versus wavenumber for 150 Hz. 
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Model 

Frequency 

(Hz) 

Eigenvalue 
Iteration 1 

Eigenvalue 
Iteration 2 

Eigenvalue 
Iteration 3 

Eigenvalue 
Iteration 4 

Background 

25 

0.095567562 

0.095840216 



Model B2 







50 

0.203297004 

0.203447193 

0.203428298 



50 

0.182848796 

0.183509603 

0.183428794 









150 



0.625676215 

0.625808477 


150 

0.617533207 

0.617683887 

0.617635310 

0.618125618 


150 

0.603655517 

0.603984475 

0.603875697 

0.604930580 


150 

0.583521008 

0.584122181 

0.583913684 








Table 3. Modal eigenvalues for the initial background geoacoustic model B1 (k,„). 


Model 

Frequency 

(Hz) 

Eigenvalue 
Iteration 1 

Eigenvalue 
Iteration 2 

Eigenvalue Eigenvalue 
Iteration 3 Iteration 4 

BtWlWCTTT ■[ 


0.096081510 

0.095852643 

0.095839143 

Model B2 






50 

0.203560099 

0.203439906 

0.203429595 


50 

0.183980197 

0.183473706 

0.183432102 i 







150 

0.625703216 

0.625677884 

0.625686586 ' 


150 

0.617737293 

0.617641628 

0.617674530 1 


150 

0.604090273 

0.603888690 

0.603957891 1 


150 

0.584284186 

0.583935380 

0.584054530 






Table 4. Modal eigenvalues for the initia 

background geoacoustic model B2 (km). 


Tables 3 and 4 list the computed eigenvalues for the initial 
background models for each iteration in the inversion process. These 
eigenvalues (k ro ), as discussed previously, are subtracted from the measured 
eigenvalues (Table 2) to form the data vector Ak(m). Generally, the lower 
frequencies (i.e., fewer modes) required fewer iterations in order for solutions 
to converge. Figures 9-11 display the results of the pertubative inversion 
technique to recover the synthetic earth sound speed profile from the initial 
background model Bl. 
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model Bl. 


As evident from the plots, the sound speed profiles are corrected 
in the proper direction relative to the background model. The profile in the 
upper tens of meters are lower in speed relative to the background, whereas 
the profile beyond about 90 m depth is virtually unaffected by the 
perturbation. Clearly, the 25 and 50 Hz signals correct the profile to greater 
depths than the 140 Hz signal. However, the 140 Hz signal approximates the 
sound speed at the water-bottom interface more accurately than the lower 
frequencies. Also apparent from the plots is that as the normal mode 
function decreases with depth, the inferred sound speed profile returns to the 
background profile. 

Figures 12-14 show the results of recovering the synthetic earth 
model using the background model B2. The results are qualitatively the 
same as the previous background model Bl. Even though the background 
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models differ significantly, the solutions converge about a sound speed profile 
which is strikingly similar to the synthetic earth model. 



Figure 12. 25 Hz Inversion of synthetic earth model from background model B2. 






Figure 14. 150 Hz Inversion of synthetic earth model from background model B2. 


The source angles in the data and the models are predominantly 
small (< 20°) and thus the perturbation technique alone can not be expected 
to resolve geoacoustic features much beyond an acoustic wavelength below 
the water-bottom interface. To resolve the deeper structures, a deep earth 
model algorithm based upon standard seismic (high source angle) signal 
processing techniques is presented in Section IV of this thesis. Since these 
methods are industry standard in the petroleum exploration industry, no 
attempt will be made to evaluate these techniques. However, the deep earth 
algorithm must be considered an essential part of the overall algorithm for 
Navy practical applications. 
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III. ACOUSTIC DATA / MODEL SUMMARY 

A. EXPERIMENTAL OVERVIEW OF PROJECT GEMINI 

A series of five shallow water, synthetic aperture array experiments 
collectively called Project GEMINI, were collected off of Corpus Christi, 
Texas in the fall of 1985. For the purpose of this thesis, only those 
measurements made at the shallow and deep sites will be discussed. The 
third site (Rubano site) has been studied extensively by previous efforts and a 
well defined geoacoustic characterization for the location is described in 
Matthews et al., 1985 and Rubano, 1980. A recent analysis of the GEMINI 
region by Duarte (1994) shows relatively good agreement between measured 
and modeled transmission loss quantities at the Rubano site. Given the 
relatively benign environment of the GEMINI region and the close proximity 
of the array sites, it was speculated by Duarte (1994) that the geoacoustic 
model, initially developed for the Rubano site, could be successfully used at 
the other two sites. However, he concluded that the TL data is particularly 
sensitive to the input geoacoustic parameters and that site-specific 
geoacoustic models are required in even in a benign, nearly range 
independent environment. 

B. ENVIRONMENTAL DESCRIPTION OF THE GEMINI AREA 

1. Sound Speed Profile 

With the exception of the shallow site, CTD and XBT casts were 
made at the commencement and completion of each tow leg. Two weeks 
prior to the experiment, Hurricane Elena swept through the Gulf of Mexico 
and, as a result, established a well-mixed water column in the GEMINI 
region. Figures 15 and 16 show the sound speed profiles for the shallow and 
deep water sites, respectively. The sound speed profile for the shallow site 
(Figure 15) indicates isothermal conditions over the entire water column. 
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Figure 

1985. 



Sound speed profile for the shallow site conducted on 11 September 


Sound Speed Profile for the Deep Site 
Sound Speed <m/») 

1335 1540 1345 1550 



Figure 16. Sound speed profile for the deep site experiment conducted on 12 
September 1985. 

The sound speed profile for the deep site (Figure 16) also shows the 
presence of a near isothermal water column down to about 41 m overlying a 
low sound speed layer of approximately 10 m in thickness. This results in a 
thermocline gradient of approximately -0.7084 s' 1 . This negative sound 
speed gradient causes sound rays to refract downward and thus accentuate 
sound propagation in the sub-bottom. 


34 



2. Bathymetry 

The Texas continental shelf is characterized by relatively smooth 
topography which gently dips toward the southeast. During Project GEMINI 
high resolution bathymetric data were collected along each of the propagation 
paths as well as at each of the station locations. Bathymetry versus range 
from the array locations are shown in Table 5. Each of the tow legs of the 
experiment was approximately 5000 m in length thus giving an average 
bottom slope at the shallow site of approximately 0.011° and about 0.054° for 
the deep site. 


Bathymetry at the Shallow 
Site 

Bathymetry at the Deep 
Site 

Range (m) 

Depth (m) 

Range 

(m) 

Depth (m) 

0 

21 

0 

62 

513 

21 

767 

61 

2288 

21 

1441 

61 

3783 

20 

2255 

60 

5460 

19 

3425 

60 



4123 

58 



4974 

58 


Table 5. Bathymetry for the GEMINI track experiments versus range. 


C. GEOACOUSTIC MODELS FOR THE GEMINI REGION 

1. Background 

In order to compare the measured TL to the SNAP and FEPE model 
estimates of TL, geoacoustic models were constructed at the shallow and deep 
water sites. A geoacoustic model describes the depth dependence of 
compressions] and shear wave speeds, compressional and shear wave 
attenuation, and density for each layer in the seafloor. In general, a 
geoacoustic model defines the true thicknesses and properties of the sediment 
and rock layers in the sea floor. Since a geoacoustic model only represents 
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the physical properties at a single point, multiple data sets must be formed to 
account for the variations in sediment properties in the horizontal. 

In a recent analysis of TL in the GEMINI area. Duarte (1994) 
showed that even in a near range-independent environment such as the 
GEMINI area, model TL fluctuations are site specific and highly dependent 
on the input geoacoustic parameters. Moreover, he found that TL data and 
model comparisons were in good agreement at the Rubano site where a site 
specific geoacoustic model was developed and incorporated into the 
propagation models. At the shallow and deep sites the Rubano geoacoustic 
data were used in the model calculations and rather poor agreement resulted. 
Thus, it appears that geoacoustic data obtained for one location can not be 
used at sites just 12 nm away, even in a seemingly range-independent 
environment. Because of this horizontal variability, it was then necessary to 
construct geoacoustic models for the shallow and deep sites of the GEMINI 
region. These are discussed in the next section. 

2. Regional Geology 

Glaciation during the Pleistocene era, 0.012-2.8 million years before 
present (mybp), created recognizable imprints on Gulf of Mexico 
sedimentation patterns. Table 6 is a generalized geologic time scale for 
reference. With the growth and subsequent retreat of the ice sheets, sea 
level fell (regressed) and rose (transgressed) accordingly. During the 
regressive sequences, numerous stream channels appeared along the 
present-day Texas continental shelf. Lowering sea levels caused rapid 
seaward progradation of river deltas, which rapidly expanded the continental 
margins. During transgressive sequences, when mean sea-level rose, vast 
areas of the former coastline were submerged. The stream channels which 
once occupied the continental shelf regions were subsequently filled with 
fine-grain sedimentary deposits. 
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Period 

Epoch 

Glacial Stage 

Sea Level 

Time 


Holocene 


Transgression 

0.012 mvbu 



Late Wisconsin 

Regression 




Mid Wisconsin 

Transgression 

0.1 mybp 



Early Wisconsin 

Regression 




Sangamon 

Transgression 




Late Illinoian 

Regression 



Pleistocene 

Mid Illinoian 

Transgression 




Yarmouth 

Transgression 




Kansan 

Regression 

2.1 mvbn 



Aftonian 

Transgression 




N'ebraskan 

Regression 

2.8 mvbn 


Pliocene 




Tertiary 




_ 52 ™y b p_ 


Table 6. Geologic time scale for the Quaternary Period (From Matthews et al., 
1985). 


The Holocene age sediment layer represents the last transgressive 
sequence in the post-glacial period. This sequence is characterized by fine 
grain silty-clays and represents the uppermost, surficial sediments in the 
GEMINI region. A thorough discussion of the sediment properties for the 
Holocene sequence is presented in Matthews et al. (1985). The thickness of 
the Holocene layer varies substantially in the GEMINI region. An isopach 
chart of the Holocene age sediments is shown in Figure 17 (Berryhill and 
Tippet, 1981). In general, the Holocene layer thickens seaward as evident by 
the seaward increase in two-way travel time from 10 to 40 ms. The dashed 
lines on the chart represent the more prominent relict stream channels along 
the shelf. The presence of stream channels on multi-channel seismic profiles 
is one means of discerning the interface between Holocene age sediments and 
the underlying late Wisconsin sedimentary sequence (Matthews et al., 1985). 
Since these stream channels were later filled by post-Wisconsin age 
sediments, the Holocene age sequences are thicker in these localized regions 
(in excess of 12 m) than in the surrounding Holocene silty-clay layers. 
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Figure 17. Sediment isopach chart of the GEMINI region expressed in 
milliseconds of two way travel time (broad solid line). Narrow solid line 
represents sediment deposited since the last low stand of sea level. Dashed lines 
represent the outline of localized thicker accumulations of sediment in ancient 
stream channels. One millisecond is approximately equal to 0.73 m (Berryhill, 
1981). GEMINI array locations are denoted by the © symbol. 


The U. S. Naval Oceanographic Office (NAVOCEANO) in conjunction 
with ARL/PSU collected 16 gravity cores in the GEMINI region (Ross et al., 
1978). The results of their analysis indicates that the Holocene age silty-clay 
layer has a mean density of 1.567 g/cm 3 and an average sound speed ratio 
(i.e., ratio of sediment sound speed to bottom water sound speed) of 0.987. 
These values differ slightly from the generic values presented in Hamilton 
(1980) for continental shelf and slope environments. However, the Ross et al. 
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data represents site-specific conditions and as such are preferred for use in 
the geoacoustic models, Unfortunately, none of the cores were able to extend 
deep enough to sample the underlying late-Wisconsin layers. Matthews et al, 
(1985) interprets this layer to be a very fine sand or sandy silt. Based upon 
the shallow water sediment tables of Hamilton (1980), the late-Wisconsin 
sediments are assumed to have a density of 1.77 g/cm ! and a sound speed 
ratio of 1.080. 

Figure 18 is a representative seismic profile (east-west transect) for 
the GEMINI region. As evident from the figure, the seismic sequence is 
characterized by gently sloping, nearly parallel reflectors broken up by 
intermittent growth faults and stream channels. The strong reflection 
pattern occurring at the Holocene-late Wisconsin interface is indicative of 
the strong impedance contrast across the interface. Since acoustic impedance 
is the product of a material’s density and sound speed, a strong reflection 



Figure 18. Multi-channel seismic reflection record across the Texas continental 
shelf near Corpus Christi. A relict stream channel from the Late Wisconsin trans¬ 
gression is shown. 


pattern on a seismic profile normally suggests a large change in lithology 
across an interface. The Holocene silty-clay layer is characterized by a 
relatively low sound speed and density whereas the late Wisconsin sandy silt 
layer is characterized by a much greater sound speed and density values. 
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This impedance contrast supports the well-defined layering structure in 
Figure 18. 

3. Geoacoustic Model Development 

Both the shallow and deep water sites have environmental and 
acoustic propagation characteristics similar to the Rubano site which was 
previously reported in Matthews et al. (1985) The sites differ only in water 
depth and in the thickness of the Holocene age sedimentary sequence. Minor 
differences in water column sound speed are also apparent, however these 
variations do not significantly impact the geoacoustic characterizations The 
shallow site is situated in about 20.5 m of water and the underlying 
Holocene-age silty clay layer is estimated to be 8 m thick. The deep site is 
situated in about 63 m of water. At the deep site the Holocene layer is 
estimated to be about 24 m thick. 

At each site the sound speed ratio for the silty-clay layer is assumed 
to be 0.987, which implies that the sound speed of the silty-clay layer is less 
that the overlying water column. Using methods described in Hamilton 
(1980), the compressional wave speed expressed as a function of depth takes 
the form of the following regression formula: 

V p = Vbw (0.987) + 1.3 s '(Z) (46) 

where Vp is the compressional wave speed in m/s, Vt>» is the bottom water 
sound speed in m/s, and Z is the depth of sedimentary layer below the 
seafloor. The product Vbw(0.987) represents the initial sound speed (V 0 ) of the 
sediment. The sediment sound speed gradient (1.3 s 1 ) was established from 
the NAVOCEANO cores and is consistent with the near-surface sound speed 
gradients reported in Hamilton (1979). Because the thickness of the 
Holocene sediments does not exceed a few tens of meters, the higher order 
terms of Equation (46) (not shown) make negligible contributions to the 
sound speed profile in the sediment layer and therefore are neglected 



(Matthews et al., 1985). As a result, the compressional wave speed in the 
Holocene layer can be expressed as a linear function of depth. 

Thermal variations at the sea floor can have a pronounced 
effect on the upper sediment layer sound speed gradient. If a periodic 
thermal fluctuation is applied at the water-sediment interface, a thermal 
wave will propagate into the sediment. The depth to which the temperature 
effects will be felt is dependent on the sediment thermal conductivity and the 
heat flow (Matthews et al., 1985). In a study of the seasonal variation of 
compressional sound speed due to temperature variability in the water 
column, Rajan et al. (1992) showed that this depth of influence was about 3 
m below the water-sediment interface. Therefore, to account for the thermal 
variability at the seafloor, Matthews et al. (1985) imposed a negative sound 
speed gradient in the upper 2.5 m of the Holocene layer. Below this depth 
the sound speed profile follows that expressed in Equation (46). 

The late-Wisconsin sediments can be modeled in a similar 
fashion. Based upon the shallow-water sediment tables of Hamilton (1980), 
these sediments are interpreted to have a sound speed ratio and density of 
1.080 and 1.77 g/cm'l respectively. In general, fine grained sands are 
characterized by large sound speed gradients near the surface decreasing 
sharply with depth (Hamilton, 1980). As a result, the sound speed profile as 
a function of depth is curve-linear. 

V p =1681 Z o° 13 (47) 

The sound speed constant in Equation (47) is based on sediment 
values at the water-sediment interface. For a buried layer, such as the late 
Wisconsin sand, compensation must be made for the lithostatic load from the 
overlying Holocene sediments. Initial values of compressional wave speed for 
the late Wisconsin layer are estimated by matching the lithostatic load 
corresponding to the overlying Holocene sediments to a similar lithostatic 
load at depth in the late Wisconsin sediment layer. At the shallow site the 
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compressional wave speed at the top of the late Wisconsin sediment layer is 
estimated to be 1717.9 m/s. For the deep site, with the increased overburden 
pressure, the compressional wave speed increases to a value of 1731.0 m/s. 
Figure 19 illustrates the sound speed profiles through the different 
lithofacies for the shallow and deep sites. Of particular interest in these 
profiles is the relative thicknesses of each propagating media. The thickness 
of the Holocene sequence more than doubles in thickness between the 
shallow and deep sites. As will be shown later, this causes significant 
fluctuations in the TL versus range curves. 
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Sound Speed Profile for the Shallow Site 


Compressions! Sound Speed (m/s) 
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Sound Speed Profile at the Deep Site 
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1500 1600 1700 1800 1900 



[ZD 


Late Wisconsin 


Figure 19. Sound speed profiles at the (a) shallow and (b) deep sites. 

For terrigenous silty-clay sediments Hamilton (1978) showed 
that the sediment bulk density varies linearly as a function of compressional 
wave speed for the first few hundred meters below the sea floor. For the 
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Holocene layer, sediment density as a function of compressional wave speed 
was computed using the equation: 

p =1.135(V P (Z)) - 0.155 (48) 

Here the regression constant (0.155) was adjusted in the equation to match 
the mean surface density of the 16 cores collected by NAVOCEANO. For the 
late Wisconsin sequence, sediment density as a function of compressional 
wave speed is expressed as: 

p=1.1175(V p )- 0.1245 (49) 

Since there were no cores which penetrated into this layer, the 
density structure is speculative and is based on historical trends from 
Hamilton (1979). 

Current estimates of low-frequency compressional wave 
attenuation have generated much debate in the acoustic community. 
Hamilton (1972) reported the results of in situ measurements of 
compressional wave speed and attenuation in various sediments off San 
Diego. These measurements and others from the literature allowed for 
analyses of the relationships between attenuation and frequency and other 
physical parameters to be made. Three important relationships resulted 
from the analyses: (1) compressional wave attenuation in marine sediments 
is approximately related to the first power of frequency in sands, muds, and 
sedimentary rocks, (2) attenuation can be estimated from sediment grain size 
(porosity) data, and (3) attenuation in marine sediments typically increases 
with depth due to the reduction of sediment porosity to a point where 
overburden pressure becomes the dominant effect. Beyond this null point, 
attenuation decreases smoothly with depth and overburden pressure 
(Hamilton, 1976). 



Using the methods proposed in Hamilton (1976), the linear 
relationship between compressional wave attenuation and frequency takes 
the form 

a = K p f " (50) 

where a is the attenuation of compressional waves in dB/m, K p is the 
attenuation in dB/m-kHz, f is the frequency in kHz, and n is the exponent of 
frequency (assumed to be 1). The proportionality constant K P at the water- 
sediment interface is estimated from the sediment porosity and mean grain 
size. According to Hamilton (1980), shallow water silty-clay sediments have 
an average porosity and mean grain size (in phi units) of 75.9% and 8.52<j>, 
respectively. These values result in an estimate of attenuation at the water- 
sediment interface which is an order of magnitude greater than published 
attenuation results of Mitchell and Focke (1980). Matthews et al. (1985) 
chose to use values for sediment porosity and mean grain size between those 
predicted for shallow-water sediments by Hamilton and those values 
indicated by Mitchell and Focke. As a result, the lower error bar estimate for 
an 8.0<j> particle was used to determine the proportionality constant. This 
yielded an attenuation value of 0.03 dB/m-kHz at the water-sediment 
interface. An expression relating the compressional wave attenuation for the 
Holocene silty-clay layer as a function of depth can thus be written as: 

K p =0.030 + 0.0016(Z) (51) 

The attenuation gradient (0.0016 dB/m 2 -kHz) was established in Matthews et 
al. (1985) and is consistent with attenuation gradients reported in Hamilton 
(1976) for silty-clays 

Compressional wave attenuation for the late Wisconsin sand 
layer, can also be empirically related to grain size and porosity. In general, 
sound attenuation for sands are much larger than those of silts, clays, and 
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muds (Hamilton, 1976). For sands, sound attenuation decreases with depth 
and increasing overburden pressure. For the late Wisconsin sequence, the 
attenuation profile is calculated from the following equation from Matthews 
et al. (1985): 


K P =0.29930 - 0.00067(Z) (52) 

Even though the Holocene and late Wisconsin age sediments 
possess enough dynamic rigidity to support shear wave propagation, shear 
wave speed and attenuation were not incorporated in the acoustic models. 
Given the shallow water depths, the relatively thin sediment cover, and the 
low frequencies involved at the GEMINI region, shear wave production was 
assumed to make a negligible contribution to the overall acoustic field (Rajan, 
1994). 

Geoacoustic model tabulations for the shallow and deep water sites 
are listed in Tables 7 and 8, respectively. As previously discussed, the 
models consist of an upper sequence of Holocene age silty-clays overlying a 
late-Wisconsin age fine sand sequence. 
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Material 

Depth (mi 

Sound Speed 

Density 

Sound Attenuation 



Vp (m/s) 

(g/ce) 

Kp (dB/m-kHz) 

Sea Surface 

0.00 

1531.93 

1 00 

0000 


1.80 

1531.96 

1.00 

0.000 


2.40 

1539.96 

100 

0.000 


3.40 

1545.48 

1 00 

0 000 


4.50 

1545.56 

1.00 

0 000 


5.50 

1545.58 

1.00 

0.000 


6.60 

1545.53 

1.00 

0 000 


7.70 

1545.55 

1.00 

0.000 


8.70 

1545.63 

1.00 

0.000 

Water 

9.40 

1545.65 

1.00 

0.000 

Column 

10.40 

1545.67 

1.00 

0.000 


11.60 

1545.68 

1.00 

0.000 


12.60 

1545.69 

1.00 

0.000 


13.20 

1545.71 

1.00 

0.000 


15.40 

1545.75 

100 

0.000 


16.70 

1545.76 

100 

0.000 


17.70 

1545.77 

1.00 

0.000 


18.60 

1545.8 

1.00 

0.000 


19.60 

1545.81 

1.00 

0 000 

Water-Sediment Interface 

20.50 

1545.82 

1.00 

0.000 


20.51 

1520.75 

1.57 

0.030 

Holocene Age 

23.00 

1515.75 

156 

0.032 

Silty-Clay 

25.50 

1523.25 

1.57 

0.030 


28.00 

1528.25 

1.58 

0.036 


28.48 

1531.13 

1.58 

0.030 


28.49 

1717.93 

1.79 

0.296 


30.50 

1721.18 

1.80 

0.293 


35.50 

1727.68 

1.81 

0.292 


40.50 

1734.18 

1.81 

0.289 


45.50 

1740.68 

1.82 

0.286 


50.50 

1747.18 

1.83 

0.284 


55.50 

1753.68 

1.84 

0.281 


60.50 

1760.18 

1.84 

0.279 


65.50 

1766.68 

1.85 

0.276 


70.50 

1773.18 

1.86 

0.274 

Late Wisconsin Age 

75.50 

1779.68 

1.86 

0.272 

Fine Sand 

80.50 

1786.18 

1.87 

0.269 


85.50 

1792.69 

1.88 

0.267 


90.50 

1799.18 

1.88 

0.265 


95.50 

1805.68 

1.89 

0.263 


100.50 

1812.18 

1.90 

0.260 


105.50 

1818.68 

1.91 

0.258 


110.50 

1825.18 

1.92 

0.256 


115.50 

1831.68 

1.92 

0.254 


120.50 

1838.18 

1.93 

0.252 


Table 7. Geoacoustic Model for the Shallow Water Site. 
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Sea Surface 0.00 

1544.74 

1.70 

1544.74 

2.50 

1544.72 

3.60 

1544.79 

4.40 

1544.85 

5.60 

1544.87 

6.50 

1544.91 

7.40 

1544.94 

8.70 

1544.98 

9.40 

1545.00 

10.70 

1545.01 

11.70 

1545.01 

12.70 

1544.98 

13.20 

1544.94 

15.60 

1544.94 

16.50 

1544.88 

17.50 

1544.85 

18.70 

1544.82 

19.60 

1544.82 

20.60 

1544.83 

21.60 

1544.85 

22.50 

1544.87 

23.60 

1544.88 

25.60 

1544.92 

26.70 

1544.94 

27.90 

1544.97 

Water 28.40 

1544.97 

Column 30.70 

1545.01 

31.70 

1545.03 

32.60 

1545.04 

35.50 

1545.09 

36.10 

1545.10 

37.40 

1545.12 

38.40 

1545.13 

39.40 

1545.15 

40.40 

1545.17 

43.70 

1544.28 

46.50 

1543.79 

47.50 

1543.66 

49.70 

1542.46 

50.50 

1541.37 

51.60 

1540.27 

52.50 

1539.44 

53.50 

1538.58 

Table 8. Geoacoustic Model for the Deep Water Site 







Material 

Depth (m) 

Sound Speed 

Density 

Sound Attenuation 



Vp (m/s) 

(g/cc) 

Kp (dB/m-UHz) 


54.60 

1538 03 

1.00 

0.000 


55.60 

1537 69 

100 

0000 


56.70 

1537.50 

1.00 

0.000 


57.60 

1537 47 

1.00 

0.000 


58.60 

1537 44 

1.00 

0.000 


59.00 

1537.40 

1.00 

0000 


60.30 

1537.40 

1.00 

0000 


61.50 

1537.31 

1.00 

0000 

Yater-Sediment Interface 

62.34 

1536.99 

1.00 

0.000 


62.35 

1516.03 

1.57 

0.030 


67.34 

1518.53 

1.57 

0.034 


69.84 

1523.53 

1.57 

0.036 


72.34 

152903 

1.58 

0.038 

Holocene Age 

7434 

1531.63 

1.58 

0.039 

Silty-Clay 

77.34 

1535.53 

1.59 

0.041 


79.84 

1538.78 

1.59 

0.043 


82.34 

1542.03 

160 

0.045 


84.84 

1545.28 

1.60 

0.046 


86.00 

1546.97 

1.60 

0.047 


86.01 

1732.07 

181 

0.289 


9100 

1738.60 

1.82 

0.287 


96.00 

1745.13 

1.83 

0.284 


101.00 

1751.66 

1.83 

0 282 


106.00 

1758.20 

1.84 

0.280 


111.00 

1764.73 

1.85 

0.277 

Late Wisconsin Age 

116.00 

1771.26 

1.85 

0.275 

Fine Sand 

121.00 

1777.79 

1.86 

0.273 


126.00 

1784.33 

1.87 

0.270 


131.00 

1790.86 

1.88 

0.268 


141.00 

1803.93 

1.89 

0.263 


151.00 

1816.99 

1.91 

0.259 


156.00 

1823.52 

1.91 

0.256 


161.00 

1830.06 

1.92 

0.254 


Table 8. Geoacoustic Model for the Deep Water Site (continued). 


D. TRANSMISSION LOSS MODELS 

In order to qualitatively assess the accuracy of the geoacoustic models 
(derived in the previous section), estimates of TL using the geoacoustic 
models as input were compared with the measured GEMINI TL data. TL 
must be considered the measure of effectiveness (MOE) for the performance 
of the inverse technique developed in this thesis. If measured and modeled 
TL agree exactly for two different geoacoustic models, then the differences in 
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the geoacoustic parameters are considered acoustically indistinguishable for 
a given frequency. 

The normal mode model (SNAP) and the finite element parabolic 
equation model (FEPE) were used in this portion of the analysis to 
quantitatively determine the impact of the perturbative technique to 
estimate of the geoacoustic model. The results of these comparisons are 
given in the next section. 

For the inverse theory problem developed in this thesis the selection of 
the proper acoustic propagation model to use is critical to the success of the 
inversion. The central issue for this study is the inversion of geoacoustic 
propertied from measured TL data assuming little to no a priori geoacoustic 
information is available. Therefore, it is crucial to have an acoustic model 
that is sensitive to even small variations in geoacoustic properties so that 
when the model and data agree well, the geoacoustic parameters are 
estimated accurately. 

Although, the GEMINI region was selected by the U. S. Navy by 
virtue of its seemingly range-independent environment, the geoacoustic 
properties were shown to be just the opposite. For this reason SNAP could 
not be used in the inverse problem because of its adiabatic assumption (i.e., 
not fully coupled). KRAKEN is similar to SNAP, except that it is a fully 
coupled normal mode model. KRAKEN was used to calculate and display the 
normal mode functions and eigenvalues (Z(m) and k ra ) that are inputs to the 
kernel of the Fredholm integral in Equation (40). The FEPE model has 
shown excellent agreement in comparison with fully coupled normal mode 
models and because FEPE is more efficient to run, it was the principle model 
used to compare inverse results with the measured TL data. 



E. TRANSMISSION LOSS MODEL/DATA COMPARISONS 

The GEMINI experiments consisted of towing a narrow-band C\V 
source at a fixed depth away from a pair of moored receivers (Lynch et al.. 
1991). An NRL J15-3 acoustic source with output tonals of 50 and 140.05 Hz 
was used during the course of the experiment. The source and receiver 
configuration is shown in Figure 20 and in Table 1. The combination of the 
two depth settings along with the two source frequencies allowed for 4 
complex pressure field magnitudes versus range measurements to be 
collected at each site, The complex pressure field measurements were 
transformed into transmission loss (TL) by the relation 
|p(r,z)| 

TL(r.z) = -20 log^-j—p (53) 

M 

where p(r,z) is the pressure field at range r and p 0 is the acoustic pressure 
field at 1 meter distance from the source. In this section the measured TL 
results from the shallow and deep sites are presented along with TL 
estimates from FEPE and SNAP. 




j 

i) 

S’u"* 

/ / !///// / // / 

Experimental source and receiver configuration for the GEMIM 


experiments. 





1. Shallow Site 

The experimental TL measurements for the shallow site are 
illustrated in Figures 21-23 along with estimates of TL based upon SNAP 
and FEPE model estimates. The 50 Hz data reveal a predominantly two¬ 
mode interference pattern whereas the 140 Hz data show a more 
complicated, multi-mode interference pattern. The measured 50 Hz TL curve 
exhibits a high degree of attenuation attributed to a low compressional 
velocity layer in the sub-bottom. The curve also exhibits a high density of 
noise spikes beyond about 2500 m in range due to local seismic profiling 
activity in the region. These rapid fluctuations should be considered as 
artifacts in the data and are not expected to be modeled accurately. 

Figures 21-23 illustrate the relatively close agreement between 
the measured and modeled TL estimates using as input the shallow site 
geoacoustic model. For the 50 Hz case (Figure 21a), beyond 2600 m there is a 
substantial increase in TL. As previously shown, the sediment thickness in 
the GEMINI region is highly variable. The sharp increase in TL is attributed 
to an abrupt increase in sediment thickness at about 2600 m in range. This 
feature is interpreted as one of the numerous late-Wisconsin age stream 
channels which criss-cross the Texas continental shelf (Matthews et al., 
1985). This type of feature would account for this significant, localized 
increase in sediment thickness. To incorporate this range-dependency in the 
TL model estimates, a second geoacoustic model was developed for the 
shallow site. Based upon the sediment isopach chart (Figure 17), the 
Holocene age silty-clay layer was increased from about 8 m to 35 m in 
thickness. As evident from Figure 23, this range-dependent geoacoustic 
model improved the 50 Hz FEPE TL estimates for ranges beyond 2600 m. 
The 140 Hz signal (Figure 21b) was not significantly impacted by this range- 
dependency since most of the 140 Hz propagation occurs in the upper few 
meters of the sediment. 
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Overall, the FEPE TL estimates are more accurate than the 
SNAP TL measurements, a result that is in stark contrast to Duarte's (1994) 
conclusions regarding TL model for which he used the Rubano site 
geoacoustic data at the shallow site. Thus, it can be concluded that the 
shallow water TL model estimates are extremely sensitive to the geoacoustic 
model used. This conclusion is further illustrated by comparing the 140 Hz 
TL model results from the shallow site (Figure 21b) with those found in 
Duarte (1994) (Figure 24). The lone difference between the two FEPE model 
estimations rests with the input geoacoustic model used. As discussed 
earlier, the shallow site geoacoustic model consists of about 8 m of near 
surface, low-velocity sediments in contrast to 17 m at the Rubano site 
(Matthews et al., 1985). The thicker low-velocity sequence (Rubano model) 
results in more loss (>10 dB) than that experienced by using the shallow site 
geoacoustic model. The SNAP TL estimate also shows reasonably close 
agreement with the measured TL data. 

Figures 21 through 23 suggest that the FEPE estimates of TL are in 
closer agreement with the measured TL field than those portrayed by SNAP. 
One possible reason for this is that SNAP is not a fully coupled normal mode 
model and, because the geoacoustic environment is range-dependent, the 
SNAP model does not permit energy to be exchanged between modes which 
the measured TL indicates is an important component. 
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<b> 


Figure 21. Comparison of experimental and theoretical transmission loss curves 
for the shallow site for a mid-column receiver. Frequency: 50 Hz (a), 140 Hz (b). 












Figure 22. Comparison of experimental and theoretical transmission loss for the 
shallow site for a bottom receiver. Frequency: 50 Hz (a). 140 Hz (b). 















Figure 23. Plot of experimental GEMINI and range-dependent FEPE model 
transmission loss for the shallow site at 50 Hz. 


MID-COLUMN RECEIVER (140 Hz) 

Range(m) 


i> 20 
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Figure 24. Comparison of GEMINI transmission loss data and FEPE transmission 
loss model estimates for the shallow site based upon the Rubano site geoacoustic 
model as input (Duarte, 1994). 
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2. Deep Site 

The deep site was occupied on two different occasions during Project 
GEMINI (10 and 12 September 1985). Because the data collected from 10 
September are incomplete only the data collected on 12 September are used 
in this study. The experimental TL measurements for the site are illustrated 
in Figures 25-26 along with model estimates of TL. As with the shallow site 
data, the measured TL is characterized by a predominant two-mode 
interference pattern for the 50 Hz signal and a multi-mode interference 
pattern for the 140 Hz signal. The FEPE and SNAP model TL estimates 
used the deep site geoacoustic model developed in the preceding section as 
input. Unlike the shallow site, the propagation path at the deep site was 
treated as range-independent. 

Although there is reasonably good agreement in amplitude between 
the modeled and measured- TL levels for the 50 Hz signal, both FEPE and 
SNAP fail to accurately model the TL fluctuations at the deep site. The TL 
fluctuations (Figure 25a) in excess of 20 dB. These fluctuations may be due 
to interfering normal modes at the site (Duarte, 1994). 

Somewhat better agreement is noted between the modeled and 
measured TL levels for the 140 Hz signal. The 140 Hz signal shows a more 
complicated multi-mode interference pattern. As illustrated in Figure 25b, 
TL fluctuations are on the order of 20 dB in TL. 
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Figure 25. Comparison of experimental and theoretical transmission loss curves 
for 140 Hz for the deep site for the mid-column receiver. Frequency: 50 Hz (a), 140 
Hz (b). 
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Figure 26. Comparison of experimental and theoretical transmission loss curves 
for 140 Hz for the deep site for the bottom receiver. Frequency: 50 Hz <a), 140 Hz 
(b). 
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IV. PERTURBATIVE INVERSION RESULTS 


Tn this section the perturbative inversion technique is applied to the 
experimental data obtained during Project GEMINI, As a result of the 
inversion, a new set of geoacoustic parameters for the shallow and deep 
water sites are estimated. As discussed previously. TL is used as the 
measure of effectiveness of the inverse geoacoustic models. Since FEPE was 
shown to be more accurate than SNAP in estimating TL, only FEPE will be 
used in the comparison between measured and modeled TL values. The 
major strength of the inversion method is that a priori information of the 
sub-bottom is not required to estimate geoacoustic parameters. In order to 
demonstrate this, we assume that the geologic and geophysical properties of 
the environment are only vaguely known in the initial background 
geoacoustic models. 

A. SHALLOW SITE 

The results of Hankel transforming the experimental pressure field 
over the entire 5000 m range at 50 and 140 Hz are shown in Figure 27. Since 
the poles of the Green's function represent the eigenvalues for discrete modes 
for the particular waveguide at the given frequency, the location of these 
poles for the bottom receiver must be coincident with pole locations for the 
mid-column receiver. In this section only the plots for the bottom receiver are 
illustrated. However, in order to form an accurate estimate of the modal 
eigenvalues of the waveguide (k me »a), the results from the two receivers are 
averaged to form a single estimate for each frequency. The eigenvalues for 
the two frequencies are tabulated in Table 9. The 50 Hz plot (Figure 27a) 
shows one distinct peak at a horizontal wave number of about 0.19 m _1 . A 
less evident, very weak pole was found to be at a horizontal wavenumber of 
about 0.17 nT 1 . For the 140 Hz case (Figure 27b), three trapped eigenvalues 
are readily identifiable from the plot. 
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s Function versus Horizontal 
Wavenumber 


0.15 0.2 0.25 

Horizontal Wavenumber (1/m) 


a 300 
■3 250 


Green's Function versus Horizontal 
Wavenumber 


0.50 0.60 0.70 0.80 

Horizontal Wavenumber (1/m) 


Frequency (Hz) 

Mode Number 

(m 1 ) 

50 

1 

0.1896861800 


2 

0.1747074050 




140 

1 

0.5541557518 


2 

0.5431438446 


3 

0.5201824925 


Table 9. Measured mode eigenvalues for the shallow site. 
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Using the modes shown in Table 9, the pertubative inverse was 
performed separately for each frequency. In order to analyze the convergence 
of the inversion, two diverse initial background models were used in the 
inversion. The first background model, BGl, consisted of a constant gradient 
(Vc = 2.0 s' 1 ) background profile with an initial sound speed of c(0) = 1500 
m/s. and with a constant density of p = 1.56 g/cm ! . The second background 
model, BG2, consisted of a constant gradient (Vc — 1.25 s' 1 ) background 
profile with an initial sound speed of c(0)=1900 m/s and with a constant 
density of p = 1.56 g/cm :i . Both background models used the experimentally 
measured sound speed profile of the water column in the inversion. The 
results of the 50 Hz inversion are displayed in Figure 28. Also shown as 
reference is the geoacoustic model developed in Section III. 


Sound Speed (m/s) 

1400 1500 1600 1700 1800 1900 





Figure 28. Results of the 50 Hz perturbative inversion using background models 
BGl and BG2. The shallow site geoacoustic model derived in Section III is also 
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Both profiles portray a lower sound speed sediment overlying a much 
higher speed geologic sequence. Even though there are similarities between 
the inverse results using BGl and BG2, it is apparent that the inverse 
results of using BG2 as the background model are in closer agreement with 
the geoacoustic model than BGl. This effect suggests that the inversion 
process is highly dependent on the initial constraints of the background 
geoacoustic model used in the inversion. When this inversion algorithm can 
be automated, numerous initial estimates will be made and the conversion 
properties can be analyzed in more detail. If only a few initial estimates are 
made, the low eigenvalues can be used to assess stability of the inversion of 
the spectral decomposition in the Backus-Gilbert solution. The higher the 
magnitude of the low level eigenvalues, the more stable and accurate the 
solution. 

Using the same initial background geoacoustic models, the results of 
recovering the sound speed profile for the 140 Hz signal are shown in Figure 
29. As the number of modes and frequency increases, large fluctuation in c(z) 
result. To remedy this situation, a simple “box-car” average filter was 
applied to c(z). As with the 50 Hz case, the 140 Hz inversion suggests a two 
layer geologic sequence consisting of a lower speed layer overlying a higher 
speed layer. 


64 




Figure 29. Results of the 140 Hz perturbative inversion using background models 
BGl and BG2. The shallow site geoacoustic model derived in Section III is also 
shown. 

It is important to emphasize that the inversion process did not 
incorporate any a priori information of the geologic and geophysical 
parameters of the seafloor. Geophysical parameters such as sediment 
density and compressional attenuation were held constant throughout the 
process. In practice however, certain knowledge of the sediment properties 
may be obtained and thus make the problem less undetermined. For 
example, by incorporating knowledge of the sound speed ratio, the profiles in 
Figures 28 and 29 can be adjusted such that the sound speed profile 
intersects the seafloor at the measured sound speed of the sediment. 
Through conventional seismic methods, information on the deep layering 
structure may be established which will place lower bounds on the 
perturbation integral. Although not used in these analyses since only low 
grazing angle data were collected, the deep earth algorithm described in a 
subsequent section could have accurately identified the sedimentary layers 
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shown in Figure 18. By incorporating this type of information into the 
inversion, the solution is less undetermined and converges rapidly to a 
solution. 

In order to test the accuracy of the shallow site inversion, model 
estimates of TL were compared with measured TL values. For the 50 Hz 
case, the results from BGl were used as the input to FEPE. It is evident 
from Figure 28 that the derived sound speed profile suggests a two layer 
sedimentary sequence consisting of a lower speed sediment overlying a 
higher speed layer. The rather sharp transition in sound speed at about 35 
m in depth indicates a rather abrupt change in lithology which we interpret 
as the change from the silty-clay layer to the sand. Therefore, generic values 
of sediment density and compressional attenuation were incorporated in the 
inverse geoacoustic model to account for this transition. For the 140 Hz case 
(Figure 29), the results from BGl were modified in a similar fashion to 
incorporate the change in density and attenuation with depth. 

Figures 30 and 31 show a comparison of the TL curves measured at 
the shallow site with the TL curves found through pertubative inversion 
technique for the same receiver combinations as in the previous section. Also 
shown in the figures are TL versus range results using the geoacoustic model 
derived in Section II. As evident from the curves, the perturbative inversion 
model results in estimates of TL that are as good or in some case better than 
using a high resolution geoacoustic model as input to the acoustic 
propagation models. 
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Figure 30. Plot of GEMINI TL data versus FEPE using as input the perturbative 
inverse solution and a site specific geoacoustic model for the mid-column receiver, 
shallow site. Frequency : 50 Hz (a), 140 Hz (b). 


67 













Figure 31. Plot of GEMINI TL data versus FEPE using as input the perturbative 
inverse solution and the site specific geoacoustic model for the bottom receiver, 
shallow site. Frequency : 50 Hz (a), 140 Hz (b). 
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B. DEEP SITE 


The Green’s function versus horizontal wavenumber plots for the deep 
site are shown in Figures 32 and 33. As with the shallow site, these plots are 
the result of Hankel transforming the experimental pressure field over the 
entire propagation path. Typically, as frequency and water depth increase, 
the number of trapped mode substantially increases. This is made apparent 
by the large number of poles of the Green’s function. The eigenvalues for the 
two frequencies are tabulated in Table 10. The 50 Hz plot (Figure 32) shows 
two distinct peaks at horizontal wavenumbers of about 0.203 rrf 1 and 0.197 
m~', A weak pole is also evident at about a horizontal wavenumber of about 
0.190 nT l . The 140 Hz plot shows four distinct modes situated at horizontal 
wavenumbers of about 0.57. 0.55, 0.54, and 0.52 m* 1 . However, only the 
latter three, more energetic modes were actually used in the inversion, since 
the pole at the horizontal wavenumber of 0.57 did not correlate with results 
from the mid-column receiver. 


Green's Function versus Horizontal 
Wavenumber 
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Figure 32. Magnitude of the depth dependent Green’s function at 50 Hz for the 
deep site experiment. Only the bottom receiver is shown. 





Table 10. Measured mode eigenvalues for the deep site. 


Using the modes shown in Table 10, the perturbative inverse was 
performed separately for each frequency at the deep site. Unlike the shallow 
site, the inversion process at the deep site was extremely sensitive to the 
initial background model used. The initial background models BG1 and BG2 
used in the previous section resulted in an unstable inverse solution. One 
explanation for this is that low order eigenvalues in the Backus-Gilbert 
decomposition result in an undetermined inverse matrix which causes the 
solution to become unstable. However, several experimental iterative runs 
showed that a background model starting with a reduced initial sound speed 
of c(0) = 1400 m/s, and a constant gradient of Vc = 2.0 s' 1 coupled with a 
constant density of p = 1.56 g/cm 3 , provided a stable sound speed profile. The 
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results of the 50 Hz inversion are displayed in Figure 34 along with the site 
specific geoacoustic model developed previously. 


Sound Speed (m/s) 

1400 1450 1500 1550 1600 1650 1700 1750 1800 



Figure 34. Results of the 50 Hz perturbative inversion at the deep site. The deep 
site geoacoustic model derived in Section III is also shown. 

Using the same initial background geoacoustic model as for 50 Hz, the 
results of the 140 Hz perturbative inversion are shown in Figure 35. In order 
to stabilize the fluctuations in c(z), a “box-car” average was once again 
applied to the 140 Hz inverse solution. Both the 50 Hz and 140 Hz sound 
speed profiles suggest that the upper sedimentary layer at the deep site 
extends deeper than indicated in the geoacoustic model. This is a valid 
result, since as discussed previously, the thickness of the Holocene age silty- 
clay layer varies substantially in the GEMINI region. 
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Figure 35. Results of the 140 Hz perturbative inversion at the deep site. The deep 
site geoacoustic model derived in Section III is also shown. 


Using values consistent with Hamilton (1978, 1980), the above sound 
speed profiles were adjusted to account for varying sediment density and 
compressional attenuation with depth. The resulting geoacoustic models 
were used as input to FEPE. Figures 36 and 37 show a comparison of the TL 
curves measured at the deep site with FEPE estimates of TL using the 
results of the corrected inversion as input. As evident from the figures, the 
estimates of TL generated by the perturbative inverse solution are extremely 
consistent with the measured TL values. 
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Figure 36. Comparison of GEMINI TL data and FEPE TL model estimates using as 
input the perturbative inverse solution and the site specific geoacoustic model for 
the deep site for a mid-column receiver. Frequency: 50 Hz (a), 140 Hz (b). 
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Figure 37. Comparison of GEMINI TL data and FEPE TL model estimates using as 
input the perturbative inverse solution and the site specific geoacoustic model for 
the deep site for a bottom receiver. Frequency: 50 Hz (a), 140 Hz (b). 
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C. DEEP EARTH SEISMIC MODEL 


The inverse technique developed in this thesis relies principally on 
shallow grazing angle (< 20° ) propagation in estimating the near surface 
geoacoustic parameters. Typically, using this inversion method, geoacoustic 
parameters can be accurately resolved to about 4/3 of an acoustic wavelength 
down from the seafloor. Inasmuch as low angle propagation paths are 
important for this inversion technique and more importantly in typical ASW 
scenarios, the technique is not effective in identifying underlying sediment 
layers beyond 4/3 acoustic wavelengths. 

To augment the capability of the inverse technique, conventional 
seismic reflection methods are proposed to resolve these deeper sediment 
structures. These methods typically rely on higher propagation grazing 
angles ( > 45 °) and thus are able to penetrate deeper into the seafloor. 
Therefore, a deep earth model is required to account for these deeper geologic 
horizons. A deep earth model was not applied to experimental data in this 
thesis; however a brief discussion is provided in this section to show how this 
information may be incorporated into the perturbative inverse technique. 

The proposed deep earth model is developed from high angle (> 45 °) 
propagation paths and low (< 25 Hz) source frequencies. The supporting data 
can be collected simultaneously as the data collected for the inverse 
technique by collecting an adequate sample of short range (i.e,, high angle) 
time-series. Subsequently, the processing then follows standard seismic 
signal processing procedures as follows: 

Step 1: The time series data are stacked or gathered using common 
midpoint (CMP) stacking algorithms. 

Step 2: The data set is deconvolved to remove artifacts from the data 
such as water borne reverberation. 

Step 3: Normal move out (NMO) corrections are applied to the CMP 
stacks. 
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Step 4: Frequency filter the data with a time dependent, bandpass 
filter. The bandpass filters are constructed in a manner to allow lower 
frequencies to pass as two way travel time increases. 

Step 5. In range-dependent sub-bottom environments, time and depth 
migration schemes should be employed to identify lateral discontinui¬ 
ties in the strata (i.e., faults, intrusions, etc.) 

Step 6: Display data in a distance versus two-way travel time plot 
such that geologic lithologies and structures can be identified. 

If steps 1 through 6 were applied to the GEMINI data, the sediment 
interface between the Holocene age silty-clay layer and the late Wisconsin 
age sand layer could be easily discerned by source frequencies as high as 50 
Hz. Incorporating the exact position of the layer into the inversion would 
reduce the undetermined state of the inverse solution and ultimately would 
improve the final product. The Navy now has the low frequency capability 
(e.g., HLF-1 sonar on submarines) on board and thus the deep earth model is 
very practical and cost effective. The deep earth model should be considered 
part of the inverse technique algorithm recommended for use in the fleet. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The technical approach in this thesis addressed two methods of estimat¬ 
ing geoacoustic parameters for acoustic propagation modeling. In the first 
part of the study propagation model estimates of TL were compared with 
measured low-frequency acoustic propagation data collected during Project 
GEMINI. Location-specific, high resolution geoacoustic models of the sub¬ 
bottom were input to the acoustic models in the first part of the study. The 
geoacoustic models derived in this section were constructed in manner 
described in Hamilton (1980) and made use of available geologic and 
geophysical data The acoustic models examined in this portion of the study 
were the SACLANTCEN normal mode acoustic propagation model (SNAP), 
and the NRL finite element parabolic equation model (FEPE). The 
comparisons were made oyer two source frequencies (50 and 140 Hz) at two 
depth settings (mid-column and near-bottom setting). Overall, it was shown 
that FEPE, using as input the detailed site specific geoacoustic 
characterizations, resulted in estimates of TL that were more accurate than 
those calculated by SNAP and were reasonably consistent with the measured 
TL field. Moreover, it was shown that TL estimates are extremely sensitive 
to the input environmental parameters. Even in a seemingly range- 
independent acoustic environment such as off the Texas coast, a small 
change in a geoacoustic parameter (e.g., sediment thickness) can result in a 
significant change in TL. 

In the second part of this study, the perturbative inversion technique 
(Rajan, 1987) was used to estimate the geoacoustic parameters for the two 
GEMINI sites. The objective of perturbative inversion is to estimate the 
depth variation in geoacoustic parameters in the sub-bottom (i.e., sound 
speed, density, attenuation, etc.) from the differences between the horizontal 
wave numbers from measured data and the corresponding wavenumbers 
computed for a background geoacoustic model. To illustrate the efficacy of 




the inverse method, the initial background geoacoustic model used in the 
inversion assumed no a priori information of the seafloor. The resultant 
inverse solution was used to predict the geoacoustic properties at each of the 
sites. As a measure of effectiveness, estimates of TL using the inverse 
results (modeled with FEPE) were directly compared to the measured 
GEMINI TL values. The final results were in excellent agreement with the 
measured TL and the resulting TL estimates were as good as or better than 
the TL estimates obtained from detailed, site-specific geoacoustic models. 

By assuming no a priori information in the inversion process and 
subsequently obtaining such excellent results in TL addresses an extremely 
critical problem for the U. S. Navy in shallow water regions. In many of the 
strategic, littoral environments knowledge of the geoacoustic properties of 
the sub-bottom is extremely limited. Direct geologic and geophysical surveys 
in these regions to determine relevant geoacoustic environmental parameters 
are both costly and time consuming and generally not consistent with the 
rapid response requirements of today’s fleet. Therefore, the use of perturba¬ 
tive inversion technique is an accurate, cost effective and practical 
alternative to obtaining acoustic performance prediction estimates in these 
regions. Data necessary for the inversion can be easily collected using 
existing Navy platforms and hardware. 

Automation of the inversion algorithm presented in this thesis could 
allow for an entire range of initial background models to be implemented in 
the inversion process. At present, one background model is assumed and 
after each iteration of the inversion, a set of eigenvalues and eigenfunctions 
are generated by the normal mode program (KRAKEN). However, in a fully 
automated mode, an entire range of background models with differing initial 
sound speeds and sound speed gradients could be analyzed by the inversion 
routine consecutively. By applying this method, a family of sound speed 
profiles can be formulated. This is especially applicable in high noise 
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environments where the eigenvalues of the inverse matrix are exceptionally 
small. 

Currently, there are numerous existing pressure field data sets (e.g., 
Bearing Stake, NATIVE I. Southwest Florida Sea Test, etc.) that the 
inversion method, as described in this thesis, can be readily applied. To test 
the robustness of the inversion, several environmental scenarios should be 
examined ranging from a weakly range-dependent environment such as the 
Southwest Florida shelf to a more geologically complex region such as the 
Northern Arabian Sea (Bearing Stake exercise). Finally, it is recommended 
that an at-sea experiment be conducted using fleet assets to collect the 
required pressure field data. 
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